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As spin glass materials have extremely slow dynamics, devious numerical methods are needed 
to study low-temperature states. A simple and fast optimization version of the classical Kasteleyn 
treatment of the Ising model is described and applied to two-dimensional Ising spin glasses. The 
algorithm combines the Pfafhan and matching approaches to directly strip droplet excitations from 
an excited state. Extended ground states in Ising spin glasses on a torus, which are optimized over 
all boundary conditions, are used to compute precise values for ground state energy densities. 



The Ising spin glass is a model for disordered magnetic 
alloys which captures the complexity of materials with 
frozen randomness and competing interactions, including 
frustration, extremely slow dynamics, and intricate mem- 
ory effects.^ The spins in the model are coupled by ran- 
dom choices of ferromagnetic or antiferromagnetic bonds, 
leading to a complex free energy landscape. There are at 
least two theoretical approaches,^ including the droplet 
and replica-symmetry-breaking pictures, used to describe 
the non-equilibrium dynamics and low-free-energy struc- 
ture of the spin glass phase space. As these theoretical 
approaches differ significantly and exact results for spin 
glasses are rare,'^ computational work has been essen- 
tial for computing scaling exponents and as a qualitative 
check of theoretical pictures.^ 

The history of the relationship between the physical 
analysis and the mathematics of numerical approaches to 
spin glasses is long and rich. In general, characterizing 
the complex free energy landscape of disordered materials 
is challenging. Direct Monte Carlo simulations are hin- 
dered by the same high free-energy barriers that inhibit 
equilibration in the physical system. It is expected^'^ 
that times t satisfying ln(t) ~ L"^ are required to equili- 
brate systems of size L, where ip > and 9 determines 
the energy scale AE{£) of excitations or domain walls on 
length scale i, AE ^ . To rephcate the many decades of 
experimental time scales and to develop a better under- 
standing of disordered systems for t ~^ oo, algorithms for 
either accelerating the approach to equilibrium or finding 
ground states in spin glasses have been developed. Many 
of these techniques (which are often generally applica- 
ble to disordered materials) are inspired by, or have in- 
spired, approaches for combinatorial optimization prob- 
lems. Parallel tempering, genetic algorithms, and ex- 
tremal optimization are examples of heuristic algorithms 
to find close approximations to equilibrated and ground 
state configurations.^ Exact general algorithms such as 
transfer matrix methods^ and branch-and-cut methods* 
require times that are exponential in powers of the system 
size, though extensive development has led to computing 
ground states in three-dimensional Ising spin glasses with 
up to 12^ spins.* 

We have found a simple algorithm for studying two- 
dimensional (2D) Ising spin glasses that combines use 
of the classical Kasteleyn city^ and application of a stan- 
dard combinatorial optimization algorithm. Besides solv- 



ing the problem on planar graphs and linking together 
these methods, we use this algorithm to study "ex- 
tended" ground states, which optimize the energy over 
choices of periodic or anti-periodic boundary conditions, 
as well as the spin configuration. This approach dramat- 
ically improves the treatment of boundary-free samples, 
so that the finite-size effects are greatly reduced. We 
have used this algorithm to determine very precisely the 
energy of the Ising spin glass in the large volume limit. 

The Edwards- Anderson Hamiltonian that is used for 
Ising spin glasses is H {{si}) ~ — JijSiSj, where the 

couplings Jij between nearest neighbor pairs of spins (ij) 
are independent identically distributed variables, fixed 
in a given sample, and the Si are Ising spin variables. 
Si = ±1, with L'^ sites i on a d-dimensional lattice. 
The distribution for Jij is generally taken either to be 
Gaussian or bimodal (the "±J" case). Barahona^*^ has 
shown that computing the ground state energy of a 3D 
spin glass (or even two coupled 2D layers) is an NP-hard 
problem.^ This implies that if the ground state of the 
3D spin glass could be efficiently computed, i.e., found 
in a time polynomial in L, many outstanding computa- 
tional problems that are believed to require worst-case 
exponential time to solve, such as the Traveling Sales- 
person Problem, could also be solved in time polynomial 
in the size of the problem. Improvements in 3D spin 
glass calculations therefore focus on reducing the numer- 
ical constants in the exponent for the expected solution 
time. 

The two-dimensional Ising spin glass (2DISG) is a case 
where exact algorithms have allowed for study of the 
ground state and density of states. These approaches 
have used two methods: the dimer-Pfaffian method (Pfaf- 
fian method) and matching to minimize frustration. 

The partition function for Ising models with arbitrary 
couplings can be solved for either open or toroidal bound- 
ary conditions by using techniques developed for the pure 
Ising model, ^ i.e., computing and summing Pfaffians, an- 
tisymmetric combinations of ordered statistical weights, 
from X sparse matrices. The ground state energy 
can be computed^^ in 0{L^) time for discrete- valued dis- 
order; the spectrum is discrete and bounded by a power 
of L. Note that the Pfaffian method uses perfect match- 
ings (dimer coverings) on a decorated lattice and requires 
a sum over four combinations of odd and even constraints 
on these matchings on a torus. 
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The fastest ground state algorithms for the 2DISG 
niap the Ising spin glass problem to the weighted perfect 
matching problem, a common problem in combinatorial 
optimization. Given a graph G = {V,E), with vertices 
V and edges E, with a weight function ui : -B — >■ M, the 
problem is to select a perfect matching, a subset of edges 
M C E where every vertex in V belongs to a single edge 
e € E, such that the total weight w{M) = J2eeM is 
minimal. The solution can be found in time polynomial 
in the number of edges. ^ Matching is the core routine 
in two mappings for finding 2DISG ground states. The 
mapping by Bieche et al}"^ uses a graph where the vertex 
set V contains the frustrated plaquettes (primitive poly- 
gons p with n^y^gpJij < 0). The edges connect points 
in V within some distance fmax- This algorithm is ex- 
act as Tmax ^ 00, but it works for a large fraction of 
cases even with small values of Tmax, especially for ±J 
disorder. Barahona's mapping^" replaces each plaquette 
with a subgraph that is connected to neighboring sub- 
graphs by dual bonds, with each dual bond crossing one 
edge in G [see Fig. 1(a)]. The subgraph edges have zero 
weight and the dual edges that cross bonds of strength 
Jij have weight Wij = \ Jij\] the subgraph comes in two 
types, assigned according to the frustration of the plaque- 
tte. These algorithms have been extremely useful, e.g., in 
studying domain walls and the nature of the ground state 
as L ^ 00.^^'-^'' Note that the graphs used are derived 
from the (sample-dependent) plaquette frustrations; this 
is not the case with our algorithm, where the graph is 
independent of the Jij , so its implementation is simpler. 

Matching algorithms have been used for planar graphs. 
The case of the torus, with periodic boundary conditions 
in both directions, has not been addressed in very large 
systems, as Pfafhan methods arc much slower (and in 
practice, mean-time exponential run- time algorithms are 
still commonly used). Studies of smaller toroidal sys- 
tems with Gaussian disorder have used the branch-and- 
cut algorithm^"^ or the transfer matrix; such studies con- 
firm that the finite-size corrections vanish much more 
quickly in toroidal geometries rather than planar geome- 
tries. It would therefore be useful to have a fast algorithm 
for finding information about the ground states for the 
2DISG on the torus. 

We have developed an approach which is not limited 
to planar graphs; it also provides significant information 
aboTit th(^ ground state on the torus. One compcment of 
this approach is a ground-state algorithm that combines 
a representation from the Pfafhan method with match- 
ing. The other component is applying this algorithm on 
the torus to find an extended ground state: the mini- 
mum energy state over all spin configurations and over 
the set of four boundary conditions (BCs). That is, we 
find the extended state ({s°} ,cr°,(j°) which minimizes 
Ti* = — Yl^ij) JijSiSjdij, with aij = 1 except on one ver- 
tical column of horizontal bonds, where = cr„, and on 
one horizontal row of vertical bonds, where aij = ah and 
ah and a^ take values ah.v = ±1. The extended ground 
state on the torus is the minimum energy state over 



the four possible combinations of BCs given by choos- 
ing (anti-)periodic BCs for each direction. The standard 
ground state for given BCs is therefore exactly found for 
J of the samples. Note that, in general, when all atj = 1, 
H* = H, so the extended ground state is equal to the 
standard ground state (this is always the case for pla- 
nar graphs, so the algorithm finds ground states of pla- 
nar graphs without modification). The extended ground 
state on the torus is also of interest in its own right. 
For example, it can be used as an edge-free background 
for studying equilibration and droplets^^ and to rapidly 
compute the energy density for large samples. 

We first give an overview of our algorithm. A spin and 
bond configuration is used to define a weighted dual lat- 
tice D which in turn is mapped to a weighted graph G. A 
minimum weight perfect matching for G is computed and 
used to identify a set of negative weight loops in D with 
the most negative total weight. These loops are exactly 
the excitations of the current configuration relative to an 
extended ground state. The configuration is thus set to 
the ground state by flipping the spins "within" each loop. 
This method can be applied to any planar graph by sup- 
plying the appropriate boundary conditions (i.e. in the 
same way as Bieche et al. and Barahona algorithms). 

A more detailed description of the method for the LxL 
toroidal square lattice starts with a list of the inputs: 
an initial conflguration c = ({sj}, a/t, cr„) and couplings 
Jij. The dual lattice D = {y,E) has edges eij e E (dj 
crosses the bond (ij) in the original lattice) connecting 
neighboring plaquettes (these make up V) on the original 
spin lattice; it also is an L x L torus. Given c, weights 
Wc for edges in E are set by Wc{eij) = JijSiSjaij; see 
Fig. 1(a). The value of the extended Hamiltonian is then 
n*{c) = -w,{E) = -Ee.,eB«'c(ei,). 

To minimize 7i*, we find the extremal (i.e., minimum 
total weight) set of negative weight loops in the dual 
graph D by computing a minimal weight perfect match- 
ing M on a related graph G. In the case of a square 
lattice, we form G by replacing each vertex in D by 
a "Kasteleyn city" subgraph, a complete graph with 4 
nodes [see Fig. 1(b)]; such mappings exist for any lat- 
tice. Weights for edges in G are zero on city edges and 
are given by Wc{eij) on edges Cij kept from D (cf. the 
Barahona algorithm, which instead uses |w(ejj)|, which 
is independent of c; frustration is incorporated via the use 
of two distinct graph decorations). Matchings in G can 
be mapped to sets of loops in D: one simply contracts 
out the Kasteleyn cities from M to arrive at a set of 
loops made of edges S C E [see Fig. l(c,d)]. The Kaste- 
leyn cities enforce the constraint that an even number of 
edges in S meet at each dual vertex (i.e. /S is a collection 
of Eulerian subgraphs of D). 

To prove the correctness of the algorithm, we first show 
that the weight of the loops that relate two configura- 
tions is proportional to the energy difference between 
the configurations. For an extended spin configuration 
c, let bij{c) = Si-Sjaij. When comparing two extended 
configurations c and c', call S the set of bonds in which 
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^ii(c) = —bij{c') (note that bij{c) = ±&y(c') always). 
Since in S, bij{c) = —bij{c') and in E\S, bij{c) — bij{c'), 
we have that 

n*{c)-n*{c') = -Y.J^Jh,{c)+Y,J^,b^Ac') 

(ij) (ij) 

= -2 ^ J^jhj{c) 

= -2 ^ w{e,,{c)), (1) 

SO that the energy difference between configurations is 
given by twice the weight of 5*. 

The proof that the minimum weight even-degree sub- 
graph always finds the ground state, then, is as follows. 
Assume, for the sake of contradiction, that there exists 
some extended spin configuration c° with a lower total 
energy than the c' returned by our algorithm from ini- 
tial configuration c. Call S the set of bonds for which 
bij{c) ~ ~bij{c'), and 5° the set of bonds for which 
b,j{c) = -hj{c°). Since H*(cO) < H*{c'), the en- 
ergy difference H*{c) - H*{c°) > n*{c) - n*{c') gives 
^SeijGSO w'c(eij) < 2X]eijes"^c(eij), which means S° 
is an even-degree subgraph of D with a more negative 
weight than S, in contradiction with the assumption that 
S is the extremal weight even degree subgraph of D. 

Note that Kasteleyn cities are often described on the 
original lattice, where loops represent a high-temperature 
expansion, but here on the dual lattice these loops con- 
tain clusters in a low-temperature expansion. The terms 
that contribute to the Pfaffian^ are products of statis- 
tical weights ±e~^'^'J over edges in loops in D and sta- 
tistical weights of unit norm from the Kasteleyn cities. 
The dominant term in the PfafHan that maximizes the 
norm of such a product minimizes the sums of the Wij 
consistent with a perfect matching in the graph G. We 
note that there has been at least one mention of using 
matching on the torus, where one of the four ground 
states was found using the Bieche et al. algorithm, but 
the utility of the extended ground state has been made 
apparent and proven by this algorithmic framework. 

This algorithm is simple to implement (given a stan- 
dard matching algorithm) and fast. On a 3.2 GHz Pen- 
tium IV processor, the extended ground state for a 100^ 
square lattice on a torus is computed in 0.8 s for Gaussian 
disorder, where we use Blossom IV'^^ for the matching 
routine. The mean solution time scales approximately as 
L^-^ through toroidal lattices of size 400^. On toroidal 
graphs with L < 128, our algorithm, which finds exact 
ground states, is at least three times as fast as our imple- 
mentation of the Bieche et al. algorithm. Note that the 
Bieche et al. algorithm does not find the exact optimal 



(a) Dual graph D (b) Dual vertex expanded to Kastelyn city 




Figure 1: (color online) An outline of the steps that convert a 
spin and bond configuration first to the dual weighted lattice 
D and then to the weighted graph G, in order to compute 
the extended ground state, (a) The original spin system, here 
with initial spins Si — 1 indicated by white circles, and bond 
configuration Jij (dashed lines) determine edge weights Wij = 
JijSiSj (taking the initial BCs to be periodic) in the dual 
graph D (solid vertices and edges), with periodic boundary 
conditions, (b) The vertices in D are replaced by Kasteleyn 
cities (light lines have zero weight in G). (c) An example set 
S of negative weight loops in D is shown, with two simple 
loops and one winding loop, (d) Heavy lines indicate the 
minimum weight perfect matching M for G (light lines are 
free edges and solid circles are vertices in G). The negative 
weight loops S (heavy dashed lines) are found by clipping out 
the Kasteleyn cities and keeping the remaining edges. Finally, 
spins are assigned by scanning across the sample: each time 
an odd number of loops is crossed, the spins are flipped (gray 
circles indicate Si — —1). In this case, the inconsistency at 
the right side is corrected by changing horizontal boundary 
conditions from periodic to antiperiodic. 



state in in all cases - in this case 1.5% of the samples 
(when r,nax — 8). Because the structure of the graph 
used in the Barahona algorithm is similar in structure to 
that of our algorithm, the two algorithms have similar 
performance, with the Barahona algorithm using slightly 
less time (about 20%) and more memory (about 20%). 

On a torus, we use this algorithm to exactly solve for 
the extended ground state, which is closely related to, 
but different than, finding the ground state of the spin 
glass for given BCs. The ground state energies for the 
four possible BCs differ by O(L^), which is the energy 
of a system-spanning domain wall. The extended ground 
state, the minimum of the four, therefore has at most an 
energy difference of 0{L^) from that for specified BCs. 

This 0{L^) difference is the same order as the expected 
finite-size correction to the ground state energy in a pe- 
riodic system, so the extended ground state is useful for 
studying energy densities. We computed the sample av- 
erage of the extended ground state energy Tig, using at 
least 5 X 10^ samples for L < 64 and at least 10^ sam- 
ples for sizes 128 > L > 64, both for Gaussian disor- 
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Figure 2: (color online) The extended ground state energy 
density eoiiJJ) for the 2D Ising spin glass on a torus is plot- 
ted vs. scaled system size L^~^ . Two scales for each disorder 
type are used, to show the linear fit at large L and the higher- 
order corrections at small L. (a,b) Assuming 6 « —0.28 for 
Gaussian disorder gives eo(oo) = —1.314788(4). (c,d) A sim- 
ilar plot using = Q for discrete values of Jij = ±1 gives 
e±(oo) = -1.401925(3). 



der (J2=l, 



Jij — 



,,y ^, 0) and for the iJ-model, Jij = ±1 

with equal probability. We then plotted the sam ple aver- 
age of the ground-state energy density, eo{L) = HqL^^, 
vs. L^~^, which will give a straight line where the lead- 
ing finite-size correction dominates. For Gaussian disor- 
der, we find a linear fit to be very good for L > 32, as 
shown in Fig. 2(a,c), for a wide range of 9 ~ —0.28(4) {9 
is not precisely determined by this method; see a sum- 
mary of results in Ref. 13) and a highly precise estimate 



eo = -1.314788(4) (cf., e.g., cq = -1.31479(2) from 
Ref. 13). Taking 9 — for the ±J data also gives a 
good fit, with = -1.401925(3) (cf. = -1.40193(2) 
from Ref. 18; finite-size effects in our L = 48 samples are 
less than those for L — 1800 samples with open BCs). 
The extra precision results from the rapid convergence 
to the thermodynamic limit in boundary- free samples, 
which can be solved much faster than standard periodic 
samples solved using branch- and-cut. 

In conclusion, we have linked together Pfaffian and 
matching methods to develop a fast algorithm for find- 
ing extended ground states in the two-dimensional Ising 
spin glass on a torus or standard ground states on pla- 
nar graphs. For many purposes, the extended ground 
states on a torus are as useful as ground states that are 
computed for a fixed choice of periodic and antiperiodic 
boundary conditions, as we show by precisely computing 
ground state energy densities. In the Pfafhan method for 
computing the partition function Z using the dual lattice 
(i.e., low temperature expansion), the dominant term in 
any of the four Pfaffians used to compute Z is due to 
this extended ground state; the partition function for a 
specified BC combination is found by carefully cancelling 
out configurations with other boundary conditions in the 
sum. Our method therefore is a combinatorial method, 
based on matching, for finding the term that dominates 
the contributing Pfaffians at low temperature. 
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